Quantum algorithm for obtaining the energy spectrum of a physical system 
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We present a polynomial-time quantum algorithm for obtaining the energy spectrum of a physical 
system, i.e. the differences between the eigenvalues of the system's Hamiltonian, provided that 
the spectrum of interest contains at most a polynomially increasing number of energy levels. A 
probe qubit is coupled to a quantum register that represents the system of interest such that the 
probe exhibits a dynamical response only when it is resonant with a transition in the system. By 
varying the probe's frequency and the system-probe coupling operator, any desired part of the 
energy spectrum can be obtained. The algorithm can also be used to deterministically prepare any 
energy eigenstate. As an example, we have simulated running the algorithm and obtained the energy 
spectrum of the water molecule. 

PACS numbers: 03.67.Ac, 03.67.Lx 



I. INTRODUCTION 

Obtaining the energy spectrum of a physical system 
is an important task in a variety of fields. In general, 
one has to solve the Schrodingcr equation of the system, 
which is a difficult task on a classical computer for large 
systems, because the dimension of the Hilbert space of 
the system increases exponentially with the size of the 
system, which is commonly defined as the number of par- 
ticles in the system. Thus, the complexity of simulating 
the quantum system grows exponentially. On a quan- 
tum computer, however, the number of qubits required 
to simulate the system increases linearly with the size of 
the system. As a result, Solving the Schrodinger equation 
of the system is more efficient on a quantum computer 
than on a classical computer fTj-Q . 

The standard quantum algorithm for obtaining the 
eigenvalues and eigenvectors of the Hamiltonian of 
a quantum system is the phase estimation algo- 
rithm (PEA) |6H13| . In the PEA, one prepares an initial 
guess state, and the algorithm randomly "selects" one of 
the energy eigenstates in the guess state and produces 
its energy as the output of the algorithm. It is worth 
mentioning here that the probability of selecting a given 
energy eigenstate is equal to the square of its overlap 
with the guess state. In reality, one is usually most in- 
terested in the energy differences between energy levels, 
instead of the absolute energy of a given energy level. In 
this paper, we present a quantum algorithm that solves 
this problem: obtaining the energy differences between 
energy levels of a quantum system described by a given 
Hamiltonian. The algorithm can also be used to prepare 
any energy eigenstate of the system. 

Our algorithm is motivated by the following observa- 
tion in simulating the dynamics of an open quantum sys- 
tem EMI: For an open system interacting with many 
environment modes, the mode that resonates with a cer- 
tain transition in the spectrum of the open system con- 
tributes the most to the decay dynamics associated with 
that transition. This property suggests a method to lo- 



cate the transition frequencies separating the different 
energy levels of a physical system. 

The basic idea of the algorithm is as follows: we cou- 
ple the quantum system to a probe qubit with a certain 
frequency, set the probe qubit in one of its energy eigen- 
states (say the excited state), evolve the whole system 
for some time, then perform a measurement on the probe 
qubit. When the frequency of the probe qubit matches 
the transition frequency between two energy levels of the 
quantum system, one observes a peak in the decay rate 
of the probe qubit. Therefore by varying the frequency 
of the probe qubit, we can locate the transition frequen- 
cies of the quantum system. We can also set the probe 
qubit to be in its ground state and measure its excitation 
dynamics. The difference is that in the former case we 
obtain the absorption spectroscopy of the system while 
in the latter case we obtain the emission spectroscopy. 

This algorithm has the following advantages: (i) There 
are several adjustable elements (initial state of the sys- 
tem, interaction operator, evolution time and system- 
probe coupling strength) that can be varied in order to 
improve the efficiency of the algorithm, (ii) The coupling 
of the system to the probe qubit can simulate a realis- 
tic interaction, and therefore the algorithm can naturally 
identify transitions that would occur in a realistic setting. 
(ill) Because of the freedom associated with choosing the 
coupling operator, the algorithm gives as an additional 
piece of output the transition matrix elements for any 
desired operator, (iv) Because the algorithm involves 
transitions between different energy eigenstates, prepar- 
ing the system in a good approximation to any particular 
energy eigenstate is less crucial than in the phase estima- 
tion algorithm. 

The structure of this work is as follows: In Sec. [n] we 
present an algorithm for obtaining the energy spectrum 
of a physical system. In Sec. IIII| we give an example 
to demonstrate the algorithm for obtaining the energy 
spectrum of the water molecule. In Sec. IIV1 we discuss 
the efficiency, the accuracy and the resource requirement 
of the algorithm, and compare our algorithm with the 
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PEA. Wc close with a conclusion section. 



II. THE ALGORITHM 

First, we make an initial guess about the range of the 
energy differences between the energy levels of the sys- 
tem, [w m i n , w max ]. We discretize this frequency range 
into j intervals, where each interval has a width of 
Aw = (w m « — U! m in) /j, and the center frequencies are 
given by uj k = u min + (k + 1/2) Aoj, k = . . . , j - 1. Wc 
now let a probe qubit couple to the quantum system, and 
we design the Hamiltonian of the whole system to be of 
the form 



H = Hs + ^wcr 2 



cAt 



(1) 



where the first term is the Hamiltonian of the system, 
the second term is the Hamiltonian of the probe qubit, 
and the third term describes the interaction between the 
system and the probe qubit. Here, ui is the frequency of 
the probe qubit (we have set h = 1), and c is the cou- 
pling strength between the probe qubit and the system, 
while o x and a z are Pauli matrices. The operator A acts 
in the state space of the system and plays the role of an 
excitation operator that transfers the initial state of the 
quantum system to another state. The frequency u> is 
taken from the frequency set w^. For a frequency of 
the probe qubit, we let the whole system evolve with the 
Hamiltonian shown in Eq. (1) for a time r. This evo- 
lution is implemented using the procedure of quantum 
simulation based on the Trotter-Suzuki formula [13]. Af- 
ter that, we read out the state of the probe qubit. Wc 
repeat the whole procedure many times in order to ob- 
tain the decay probability. Then we change the probe 
frequency and repeat this procedure until we cover all 
the frequencies in the range [w m ; n , w max ]. 

Setting the probe qubit in its excited (ground) state, 
when the frequency of the probe qubit matches the tran- 
sition frequency between two energy levels of the quan- 
tum system, the probe qubit has the fastest decay (exci- 
tation). For example, in the case where the initial state 
of the probe qubit is the excited state, the final state of 
the probe qubit is: 

Pp {r) = Tr s [[/(r) (|^)(V S | ® |1>(1|) U\t)], (2) 

where Trg]---] means tracing out the system degrees 
of freedom. The unitary evolution operator U(t) = 
cxp (— iHt), H is given in Eq. (1), \tp s ) is the initial state 
of the system, and |1) represents the excited state of the 
probe qubit, while |0) represents the ground state of the 
probe qubit. The quantity of interest to us now is the 
decay probability of the probe qubit Pdecay = (0|jO p (r)|0). 
By plotting Pdecay as a function of the probe-qubit fre- 
quency, we can obtain the absorption spectrum of the 
system. If there arc no degeneracies in the transition fre- 
quencies, at most one transition (denoted by i — > j) in 
the system will contribute to the decay dynamics of the 
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FIG. 1: Quantum circuit for obtaining the energy spectrum of 
a physical system. The first input register represents a probe 
qubit, and the second input register represents the system 
whose spectrum we are trying to obtain. 



probe qubit (taking into consideration the possibility of 
degenerate transitions makes the derivations longer but 
does not affect our main results). In this case, we obtain 
the result 



P, 



decay 



(Ej -Ei- w fc ) 



where Qij 



2c|(</5,|,% i )|, and Oi 



(3) 



Qij + (Ej - Ei - uj k ) . is the i-th energy 

cigenstate of the system and Ei is the corresponding 
eigenenergy. Eq. (3) describes Rabi-oscillation dynamics, 
where the system and probe exchange an excitation. The 
second factor on the right-hand side is the maximum 
oscillation probability, and it depends on the relation 
between the matrix element for a given transition and 
the system-probe detuning for that transition. The 
third factor is the overlap between the initial state of 
the system and a given energy eigenstate. 

In general, the interaction between the probe qubit and 
the quantum system should be weak such that the widths 
of the peaks are small and one obtains accurate results. 
The evolution time r should ideally be large (ct ~ 1), 
such that the change of the system is clear and the peaks 
in the spectroscopy have high resolution. 

The procedure of the algorithm is as follows: (i) pre- 
pare a quantum register i?s, which encodes the state of 
the system, in state \ip s ), and the probe qubit in state |1); 
(ii) implement the unitary operator U (r) = exp (—iHr) 
where H is given in Eq. (1); (Hi) read out the state of 
the probe qubit; (iv) perform steps (i) - (Hi) many times 
in order to obtain good statistics and calculate the decay 
probability; (v) repeat steps (i) - (iv) for different fre- 
quencies of the probe qubit. From the above procedure, 
one obtains the absorption spectroscopy of the system. 
One can also set the probe qubit in its ground state |0), 
and perform the above steps in order to obtain the emis- 
sion spectroscopy. The quantum circuit for steps (i) - 
(Hi) is shown in Fig. 1. 
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III. EXAMPLE: OBTAINING THE ENERGY 
SPECTRUM OF THE WATER MOLECULE 

In the following, we present an example that demon- 
strates how the algorithm would perform in obtaining the 
energy spectrum of the water molecule. 

To apply the quantum algorithm presented above, first 
we have to map the state space of the water molecule to 
the state space of the qubits. Using the mapping tech- 
nique introduced in Ref. considering the C 2 v an d 
1 Ai symmetries known from quantum chemistry of the 
water molecule, we can minimize the number of qubits 
needed to represent the water molecule on a quantum 
register. Note that the symmetries can be used with 
the PEA in order to optimize the algorithm in the same 
way that we have used them in the water-molecule ex- 
ample. These symmetries would not lead to an expo- 
nential speedup in cither algorithm. In other words, 
knowing the symmetries is not crucial for running the 
algorithm and for having a polynomial scaling of re- 
sources. The Hamiltonian for the water molecule is given 
in Ref. [lH and shown below. The Hartree-Fock wave 
function for the ground state of the water molecule is 
(l ai ) 2 (2ai) 2 (16 2 ) 2 (3ai) 2 (16i) 2 . Using the STO-3G ba- 
sis set [l8| and freezing the first two a\ orbitals, we 
construct a model space with 1 A\ symmetry that in- 
cludes the 3ai,4ai, I61 and lb 2 orbitals and we consider 
only single and double excitations to the external space 
for performing the multi-reference-configuration interac- 
tion (MRCI) calculation. The MRCI space is composed 
of 18 configuration state functions. Therefore at least 5 
qubits are required to represent the state of the water 
molecule in this calculation. 

In order to optimize the implementation of the algo- 
rithm, it is useful to have a priori knowledge of the molec- 
ular states and their symmetries. This can be done using 
quantum-chemistry algorithms on a classical computer. 

The Hamiltonian of the water molecule in the form of 
second quantization is 

H = J2(p\T + V N \q)ala q ~ £ (p |(g \V e \ r)\ s) 

p,q p,q.r,s 

(4) 

where \p) is the one-particle state, a|, is its fermionic cre- 
ation operator, and T, Vn, and V e are the one-particle 
kinetic operator, nuclear attraction operator and the two- 
particle electron repulsion operator, respectively. 

For the initial state, we prepare the system register 
in the simple state |00010), which is close to the true 
ground state. Then we implement the unitary operation 
U = exp (—iHr). For the interaction operator A, we set 

A = (A 1 +A 2 + A 3 + A 4 + A 5 )/V5, (5) 

where A\ ^I®I®I®I®a x , A 2 = I <g> I ® I ' ® a x ® I, 
A 3 = I®I®a x ®I®I, A 4 = I(g>a x (g>I(g)I(g)I, 
A5 =(T x (g)I(g)I(£)I<g)I. We set the coupling strength 



c = 0.005 and the evolution time r = 500 (here we mea- 
sure energies in units of Hartree and time in units of 
Hartree -1 ). We vary the frequency of the probe qubit 
in the range to £ [0.4, 2.0], which is divided into 200 in- 
tervals, and run the circuit shown in Fig. 1. We obtain 
the spectrum shown in Fig. 2 for the transition frequen- 
cies between the ground state and several excited states. 
From the figure we can see that the spectroscopy obtained 
using our algorithm is in very good agreement with the 
known transition frequency spectrum (in red) of the wa- 
ter molecule. 

The coupling strength c and the evolution time r can 
be adjusted to improve the resolution of the peaks and 
the accuracy of the results. In order to demonstrate this 
point, we now set c = 0.001 and r = 2500. We focus on 
the second and the third peaks as shown in the inset of 
Fig. 2. We can see that the widths of the peaks are re- 
duced and the resolution of the peaks is now higher. We 
also observe a small peak at the frequency of the transi- 
tion between the second and the eighth energy levels. 

From Fig. 2, we can see that some transitions between 
the ground state and the excited states are barely visible. 
Their decay probabilities can be improved by construct- 
ing a different operator A. The choice for the operator 
A in Eq.(5) includes single-qubit operators with all the 
qubits represented. With this choice most of the desired 
resonance peaks are observed in the simulation. However, 
as can be seen in Fig. 2, some peaks are very low. We use 
two-qubit operators in order to look for any such "miss- 
ing" peaks. We have tried a few different choices, and 
we only show the one that resulted in all the peaks being 
visible. In principle, even if no single operator (as hap- 
pened in our example) produces all the resonance peaks, 
one can still construct the spectrum by putting together 
the information obtained from the different choices for 
A. We define the operators A§ = I ®1 ®I ®a x ®a x and 
At = I ® I ® cr x ® a x ® /, set the interaction operator 

A = (A 1 + A 2 + A 3 + A 6 + A 7 )/V5, (6) 

and run the algorithm. The results are shown in Fig. 3. 
We can see that now all the expected peaks are clearly 
s visible. 

In our algorithm, we can transfer the initial state of the 
system to another state through the interaction operator 
A. Therefore the initial state of the system is not of 
crucial importance to the success of the algorithm. Here 
we give an example that demonstrates how the PEA can 
fail when the initial state is not a good approximation of 
the desired state, but where our algorithm still succeeds. 

In the PEA, the success probability of the algorithm 
depends on the overlap of the initial guess state with the 
desired eigenstate of the system. In the previous exam- 
ple, if the initial state of the water molecule is prepared 
in state 1 11111), the overlap of this state with any of the 
18 eigenstates (in our example, the dimension of the state 
space of the water molecule is 18) of the water molecule 
is zero. Therefore the PEA will fail in such a case. Our 
algorithm, however, still works. 
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FIG. 2: (Color online) Transition frequency spectrum be- 
tween the ground state, \(p Q ), and the first ten excited states 
(l^j), i = 1,2, ... , 10) of the water molecule. The blue solid 
curve represents the decay probability of the probe qubit at 
different frequencies with the coupling coefficient in Eq. (1) 
c = 0.005 and the evolution time r = 500, and the operator 
A as shown in Eq. (5). The red dotted vertical lines represent 
the known transition frequencies between the ground state 
and the first ten excited states of the water molecule. In the 
inset, the second and the third transition frequencies shown in 
blue were obtained using c = 10~ 3 and r = 2500. The green 
vertical dashed line represents the known transition frequency 
(1-7) between the first and the seventh excited states. 



0.9 



o 
o 

!5 

CD 

o 

Q_ 

>, o.: 

co 

o 

CD 

Q 



0.0 a-** 

0.4 




0.8 1.2 1.6 

Frequency a> (hartrees) 



2.0 



FIG. 3: (Color online) Same as in Fig. 2, except the operator 
A is set as shown in Eq. (6). 



We set the operator A to be 
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(7) 



where B\ = a x ®a x ® cr x 
a x , B 3 = a x <g)a x <E)I<E>a :t 
B5 = a x eg a x (g a x eg I < 
B 7 = a x eg I eg a x ® I £ 
B() = a x ® o x ® a x ® ( 



^a x ®I, B 2 = cr x ®a x ®a x ®I® 
g)(7 x , i?4 = <7 x ®I®a x ®<7 x ®a x , 
) 7, _B 6 = a x eg 7 eg I eg <j x eg cr x , 
/, B 8 = a x ® I ® I ® I ® a x , 
x eg 0%. And we use the state 
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FIG. 4: (Color online) Transition frequency spectrum be- 
tween the first 11 eigenstates of the water molecule and the 
state 1 11111). The blue solid curve represents the decay of 
the probe qubit at different frequencies with the coupling co- 
efficient shown in Eq. (1) c = 0.002 and the evolution time 
r = 800 and operator A as shown in Eq. (7), in simulating 
the algorithm, and the red dotted vertical line represents the 
known eigenenergies of the first 11 eigenstates of the water 
molecule. 



1 11111) as the initial state of the system. We set the 
coupling coefficient c = 0.002, and the evolution time 
r = 800, and run the algorithm. The results arc shown 
in Fig. 4. From this figure we can see that the algorithm 
still has a high success probability in obtaining the energy 
spectrum of the water molecule. 



IV. DISCUSSION 

In the following, we discuss the factors that affect the 
efficiency, the accuracy and the resource requirements of 
the algorithm, and we compare our algorithm with the 
phase estimation algorithm. 

The efficiency of the algorithm is most naturally de- 
fined through the number of times that the circuit in 
Fig. 1 must be run in order to identify the peaks in the 
spectrum. This number is proportional to the number of 
frequency points that need to be used and the number 
of times that the circuit needs to be run for a single fre- 
quency. Most physical systems have typical energy scales 
that are linear in the system size (for the total energy), 
while some unusual systems exhibit a polynomial depen- 
dence with relatively small exponents. The energy scale 
for the low-energy spectrum might be even smaller than 
that scale. The number of frequency points that need 
to be used in the algorithm, which is proportional to 
the frequency range, therefore scales polynomially with 
the system size. The number of times that the circuit 
needs to be run for a single frequency must be at least 
proportional to 1/Pdecay in order to observe a peak. It 
should also be mentioned here that each single run of 
the algorithm is essentially a quantum simulation of the 
dynamics, which scales polynomially with the size of the 
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system [l9j . 

We note here that, for large systems, there is an ex- 
ponentially large number of energy eigenstates, and de- 
termining the entire spectrum of a large system exhibits 
exponential complexity. However, one is usually not in- 
terested in all of the energy eigenstates, but rather a very 
small fraction of them, that is, a polynomial number of 
energy eigenstates. The energy eigenstates of interest 
could for example be the low-lying energy levels or the 
energy levels that are connected with the ground state 
by strong electric-dipole transitions. Once the criterion 
for the "energy levels of interest" is specified, and their 
number is small (or at least not exponentially large), the 
complexity of the algorithm does not grow exponentially 
with the size of the system any more. The part of the 
spectrum of interest will appear naturally in our algo- 
rithm, because the system will undergo transitions that 
mimic those of the simulated system. 

Since the heights of the peaks depend on the product 
QijT, we assume as a "worst-case scenario" that Qyr <C 
1 and find that the decay probability in Eq. (3) at the 
center of a given peak (i.e. at cjj. = Ej — E{) can be 
approximated as 

Pdcc ay «c 2 r 2 |(^.|A|^)| 2 |(^|^)| 2 . (8) 

From the above equation, we can see that the decay prob- 
ability, and therefore the efficiency of the algorithm, de- 
pends on the coupling strength, the evolution time, the 
interaction operator and the initial state of the system. 
Note that we can also use a number of qubits in paral- 
lel as probe qubits in order to improve the efficiency of 
the algorithm. We now address the roles of the different 
factors appearing in Eq. (6). 

The parameters c and r define the accuracy of the algo- 
rithm. The accuracy is given by the width of the peaks, 
and this width is given by max [c(y>j|A|y>j), 1/r] [15j. To 
obtain accurate results, we need to set c to be small so 
that we have weak system-probe coupling and the evo- 
lution time t is set to be large such that the change of 
the system remains observable. Note that the accuracy 
is set by the experimenter, independently of the size of 
the system. It is also worth noting here that the size of 
the frequency intervals Au is set by the choice of c and 
t: Aoj should be smaller than the width of the peaks in 
order to avoid missing some of the peaks, but there is no 
point in reducing Aw far beyond this point. 

Since Pdccay depends on (tp^\A\tp^j, the algorithm can 
also be used to obtain the matrix element (ip-\A\ip^j for 
any operator A and any two energy eigenstates, provided 
that this matrix clement is not exponentially small. For 
this purpose, it would be ideal to set the initial state to 
one of the states j^J, which can be achieved as will be ex- 
plained shortly. One can then use the height of the peak 
to obtain {tp j\A\p .j) . In this context it is also worth noting 
that since the height of the peak depends on (tp^A^^), 
certain transitions might not result in any peaks if the 
relevant matrix element vanishes. By designing A to be 
a physically-relevant operator, e.g. the operator that de- 



scribes the coupling of a molecule to an electric field, one 
can identify transitions that would occur under electro- 
magnetic irradiation of a molecule. Needless to say, A is 
not restricted to be a naturally occurring operator. 

The last factor in Eq. (6) is the overlap between the ini- 
tial state and any given energy eigenstate (which serves 
as the initial state in a given transition). In principle, 
preparing an initial state that has a large overlap with 
any given energy eigenstate can be a difficult task, possi- 
bly involving exponential scaling in the size of the system. 
However, a crucial point here is that once we observe a 
transition at the end of given run of the algorithm, we 
know that the final state of the system is the final en- 
ergy eigenstate of the relevant transition. We can now 
use this state as the new initial state and rerun the al- 
gorithm. If the new initial state is different from a state 
that we wish to examine, we can convert the state of the 
system into the desired state by adding or subtracting 
from the system the required energy difference, which 
we would know at least approximately [20j . Even if the 
two states in question were macroscopically different, it 
should be possible to go from one of them to the other 
via an at-most polynomially large number of transitions 
each of which involves single- or few-body operators. In 
many cases of practical interest, a relatively small num- 
ber of transitions are needed to connect the energy levels 
of interest, e.g. the low-lying energy eigenstates. 

We note here that, since our algorithm relics on 
changes in the energy to keep track of the state of the 
system, one cannot tell whether a given energy level is de- 
generate or not, and in the case of degeneracy one cannot 
tell which final state is obtained upon detecting the rel- 
evant transition. If one wishes to check for degeneracies, 
one could add a few small perturbations to the Hamilto- 
nian of the system, and for most physical systems these 
perturbations will lift the degeneracies in the spectrum. 

Finally, we compare our algorithm with the PEA. In 
the PEA, one prepares an initial state that ideally has 
a large overlap with the desired energy eigenstate, and 
the algorithm produces the energy of that state. In cases 
where the desired energy eigenstate has a complicated 
form or whose form is unknown, it can become impossible 
to prepare a guess state that has any substantial overlap 
with the desired state. The algorithm would fail in this 
case. In our algorithm, the initial state does not need 
to have a large overlap with any particular state. As 
mentioned above, the observation of a transition in the 
probe signals a corresponding transition in the system. 
The post-transition state can now be treated as the initial 
state for the next step in the algorithm. This way, one 
can guide the system to any energy eigenstate, including 
the ground state. The freedom in choosing the operator 
A allows additional controllability for this purpose. 

We note here that there is no single choice of the oper- 
ator A that is needed in order to obtain a certain energy 
difference, say between the ground state and first-excited 
state. As explained above, it should be possible to go 
from any state to any other state via a relatively small 
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number of intermediate states. An exception might be 
glassy and similar frustrated systems with a large number 
of vastly different low-energy states. However, there is no 
known efficient algorithm, classical or quantum, for ex- 
haustively identifying the low-energy states of such com- 
plex systems. 

In the PEA, one needs to have a good idea about the 
form of the energy eigenstates of interest. In our algo- 
rithm no such a priori knowledge is needed. If one is in- 
terested in the low-energy spectrum, the relevant states 
would show up naturally in the spectrum. This prop- 
erty is demonstrated in the example presented in Sec. IIIIl 
showing that our algorithm can work in cases where the 
PEA fails. 

In the PEA, one obtains the absolute eigenenergy of 
the system. For a large system, the absolute eigenenergy 
could be a large number, much larger than the separation 
between the energy eigenstates of interest. This large 
overall energy would appear as part of the output, thus 
taking up resources such as additional index qubits. In 
our algorithm, one obtains the energy difference between 
energy levels, therefore avoiding the unnecessary readout 
of any overall energy shift. Note that the number of 
qubits required for implementing our algorithm is the 
same as that in the optimized version of the PEA [H, [2l[ . 



V. CONCLUSION 

We have presented a hybrid analogue/digital quantum 
algorithm for obtaining the energy spectrum of a physical 
system. The algorithm provides more flexibility than the 
phase estimation algorithm. It can also be used to sim- 
ulate a realistic interaction, and naturally identify tran- 
sitions that would occur in a realistic setting. The al- 
gorithm can also be used to prepare any desired energy 
eigenstate of a physical system. 
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